A  minimum  entropy  principle  of  high  order  schemes  for  gas  dynamics 

equations1 

Xiangxiong  Zhang2  and  Chi- Wang  Shu3 

Abstract 

The  entropy  solutions  of  the  compressible  Euler  equations  satisfy  a  minimum  principle 
for  the  specific  entropy  [10].  First  order  schemes  such  as  Godunov-type  and  Lax-Friedrichs 
schemes  and  the  second  order  kinetic  schemes  [6]  also  satisfy  a  discrete  minimum  entropy 
principle.  In  this  paper,  we  show  an  extension  of  the  positivity-preserving  high  order  schemes 
for  the  compressible  Euler  equations  in  [13,  14],  to  enforce  the  minimum  entropy  principle 
for  high  order  finite  volume  and  discontinuous  Galerkin  (DG)  schemes. 


AMS  subject  classification:  65M60,  76N15 

Keywords:  compressible  Euler  equations;  discontinuous  Galerkin  method;  high  order  ac¬ 
curacy;  gas  dynamics;  finite  volume  scheme;  essentially  non-oscillatory  scheme;  weighted 
essentially  non-oscillatory  scheme;  minimum  entropy  principle 

1Research  supported  by  AFOSR  grant  FA9550-09-1-0126  and  NSF  grant  DMS-1112700. 

2Division  of  Applied  Mathematics,  Brown  University,  Providence,  RI  02912.  E-mail:  zhangxx@dam.brow- 
n.edu 

3Division  of  Applied  Mathematics,  Brown  University,  Providence,  RI  02912.  E-mail:  shu@dam.brown.edu 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

18  JUL  2011 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

A  minimum  entropy  principle  of  high  order  schemes  for  gas  dynamics 
equations 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Brown  University, Division  of  Applied  Mathematics, Providence, RI, 02912 


3.  DATES  COVERED 

00-00-2011  to  00-00-2011 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

TR-2011-19 


9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

The  entropy  solutions  of  the  compressible  Euler  equations  satisfy  a  minimum  principle  for  the  specific 
entropy  [10].  First  order  schemes  such  as  Godunov-type  and  Lax-Friedrichs  schemes  and  the  second  order 
kinetic  schemes  [6]  also  satisfy  a  discrete  minimum  entropy  principle.  In  this  paper,  we  show  an  extension 
of  the  positivity-preserving  high  order  schemes  for  the  compressible  Euler  equations  in  [13, 14],  to  enforce 
the  minimum  entropy  principle  for  high  order  finite  volume  and  discontinuous  Galerkin  (DG)  schemes. 

15.  SUBIECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

20 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1  Introduction 


The  one  dimensional  version  of  the  compressible  Euler  equations  for  the  perfect  gas  in  gas 
dynamics  is  given  by 


with 


wt  +  f(w)x  =  0,  t  >  0,  x  G  R, 


m 


pu~  +  p 
(E  +  p)u 


m  =  pu, 


E 


2  <M  + 


V 

7-1’ 


(1.1) 


where  p  is  the  density,  u  is  the  velocity,  m  is  the  momentum,  E  is  the  total  energy  and  p 
is  the  pressure.  We  consider  the  initial  value  problem  for  system  (1.1)  with  the  initial  data 


w0  (a). 


It  is  well  known  that  entropy  inequalities  should  be  considered  for  general  hyperbolic 
conservation  laws.  The  generalized  entropy  function  for  (1.1)  is  a  smooth  convex  function 
U( w)  with  an  entropy  flux  F( w)  such  that  the  following  relation  holds: 


U1  f  —  F 

WWAW  ±  w  • 


Entropy  solutions  of  (1.1)  are  weak  solutions  which  in  addition  satisfy  U(w)t  +  F( w  )*  <  0 
in  the  sense  of  distributions  for  all  the  entropy  pairs  ( U ,  F). 

In  [10],  a  minimum  principle  of  the  specific  entropy  S(x,t)  =  In  was  proved  for  the 
entropy  solutions: 

S(x,  t  +  h)  >  min{S'(|/,  t)  :  \y  —  x\  <  \  |w||oo/i}. 

The  first  order  schemes  including  Godunov  and  Lax- Friedrichs  schemes  preserve  a  similar 
discrete  property  [10].  In  [6],  a  first  order  kinetic  scheme  for  multi- dimensional  cases  on  a 
general  mesh  and  a  second  order  kinetic  scheme  satisfying  the  same  property  were  discussed. 
However,  it  seems  difficult  to  construct  higher  order  minimum-entropy-principle-satisfying 
schemes. 
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In  this  paper,  we  will  discuss  the  minimum  entropy  principle  of  an  arbitrarily  high  order 
scheme  on  a  rectangular  or  unstructured  triangular  mesh.  To  have  the  specific  entropy 
well-defined,  the  very  first  step  is  to  guarantee  the  positivity  of  density  and  pressure  of  the 
numerical  solution,  which  can  be  done  for  a  high  order  finite  volume  or  a  discontinuous 
Galerkin  (DG)  scheme  following  [7,  13,  14,  15].  The  main  idea  of  positivity-preserving 
techniques  for  high  order  schemes  in  [13]  is  to  find  a  sufficient  condition  to  preserve  the 
positivity  of  the  cell  averages  by  repeated  convex  combinations,  namely, 

1.  Use  strong  stability  preserving  (SSP)  high  order  time  discretizations  which  are  convex 
combinations  of  Euler  forward.  For  more  details,  see  [9,  8,  4,  3].  Then  it  suffices  to 
End  a  way  to  preserve  the  positivity  for  the  Euler  forward  time  discretization  since  the 
set  of  states  with  positive  density  and  positive  pressure  is  convex. 

2.  Use  first  order  schemes  which  can  keep  the  positivity  of  density  and  pressure  as  building 
blocks.  High  order  spatial  discretization  with  Euler  forward  is  equivalent  to  a  convex 
combination  of  formal  first  order  schemes,  thus  will  keep  the  positivity  provided  a 
certain  sufficient  condition  is  satisfied. 


3.  A  simple  conservative  limiter  can  enforce  the  sufficient  condition  without  destroying 
accuracy  for  smooth  solutions. 


In  fact,  the  methodology  above  can  be  used  to  enforce  any  property  for  high  order  schemes 
as  long  as  the  states  satisfying  this  property  form  a  convex  set.  In  particular,  we  will  show  in 
Section  2  that  the  specific  entropy  function  is  quasi- concave,  thus  the  following  set  is  convex, 


G  =  <  w  =  m  p  >  0,  p  >  0,  and  S  >  S0  =  min  S^wotT)) 

1  \  E  ) 


Therefore,  we  can  easily  derive  a  sufficient  condition  for  a  high  order  scheme  to  keep  numerical 


solutions  lie  in  G,  i.e.,  the  minimum  of  the  specific  entropy  at  any  later  time  will  be  bounded 
from  below  by  the  initial  minimum.  Then  a  straightforward  extension  of  the  limiter  in  [13] 
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can  enforce  this  condition  without  destroying  conservation.  This  limiter  will  not  destroy 
accuracy  for  generic  smooth  solutions,  to  be  explained  in  Section  2. 

The  conclusion  of  this  paper  is,  by  adding  a  simple  limiter  which  will  be  specified  later 
to  a  high  order  accurate  finite  volume  scheme,  i.e.,  the  essentially  non-oscillatory  (ENO) 
and  the  weighted  ENO  (WENO)  finite  volume  schemes,  or  a  discontinuous  Galerkin  scheme 
solving  one  or  multi-dimensional  Euler  equations,  with  the  time  evolution  by  a  SSP  Runge- 
Kutta  or  multi-step  method,  the  final  scheme  satisfies  the  minimum  entropy  principle  and 
remains  high  order  accurate  for  generic  smooth  solutions. 

The  paper  is  organized  as  follows.  We  first  describe  the  one-dimensional  case  in  Section  2. 
Then  we  discuss  the  two-dimensional  cases  in  Section  3.  In  Section  4,  we  show  the  numerical 
tests  for  high  order  DG  schemes.  Concluding  remarks  are  given  in  Section  5. 

2  The  one-dimensional  case 

2.1  Preliminaries 

Lemma  2.1.  S'(w)  =  In  is  a  quasi-concave  function,  namely,  the  following  inequality 
holds, 

S( AiWi  +  A2w2)  >  min{S’(wi),  5(w2)},  if  pi,  p2  >  0, 
where  wq  7^  w2  Ai,  A2  >  0  and  X±  +  A2  =  1. 

Proof.  Let  U( w)  =  —ph(S( w)),  then  I/ww  is  positive  definite  if  and  only  if  p(h'(S)  — 
7 h''(S))  >  0  and  h'(S)  >  0,  see  [5].  In  particular,  we  can  take  h(S)  =  S.  Let  w  = 
A1W1  +  A2w2  and  S*  =  min{S'(wi),  S(w2)}.  I/ww  >  0  implies  U( w)  <  Aif/(wi)  +  A2G(w2). 
So 

-pS( w)  <  -AipiS'(w1)  -  A2p2S'(w2)  <  -X1p1S*  -  X2P2S*  =  - pS *. 

Thus  we  have  S( w)  >  S*  =  min{S'(wi),  5(w2)}.  □ 

Lemma  2.2.  For  a  vector  valued  function  w(x)  =  (p(x),  m(x),  E(x))T  defined  on  an  interval 
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Ij  =  [Xj_i,Xj+i\  satisfying  p(x)  >  0  for  all  x  G  Ij,  we  have 


S 


1 

Ax 


w (x)dx  )  >  min>S'(w(x)), 


is/,- 


where  Ax  =  x,- ,  1  —  x,-_i. 

J~r2  J  2 


Proof.  Define  U (w)  =  —  pS  then  U  is  a  convex  function.  Let  p  =  A-  p(x)dx  and  w  = 
A-;  fj  w (x)dx.  Jensen’s  inequality  implies 

\ 

pS( w)  =  U  I  — —  /  w (x)dx  I  <  — [  U  (w (x))dx 


Ax 


Ax 


=  —  — — p(x)S(w(x))dx  <  —  pmmS(w(x)). 


Ax 


x£lj 


□ 


Therefore,  G  defined  in  (1.2)  is  a  convex  set  by  the  concavity  of  pressure  and  Lemma 
2.1.  The  entropy  solutions  are  in  the  set  G,  see  [10]. 

Consider  a  first  order  scheme  for  (1.1) 


W"+1  =  Wn  -  A[f(wn,  w"+1)  -  f(w”_!,  w 


(2.1) 


where  f(-,  •)  is  a  numerical  flux,  n  refers  to  the  time  step  and  j  to  the  spatial  cell  (we  assume 
uniform  mesh  size  only  for  simplicity),  and  A  =  is  the  ratio  of  time  and  space  mesh 
sizes,  w f  is  the  approximation  to  the  cell  average  of  the  exact  solution  v(x,  t)  in  the  cell 
Ij  =  [Xj_i,Xj+i],  or  the  point  value  of  the  exact  solution  v(x,t)  at  Xj,  at  time  level  n.  For 
Godunov,  Lax-Friedrichs  and  kinetic  type  fluxes  [6],  the  scheme  (2.1)  satisfies  that  w"  being 
in  the  set  G  for  all  j  implies  the  solution  w ”+1  being  also  in  the  set  G.  This  is  usually 
achieved  under  a  standard  CFL  condition 


A  ||  (M  +  c)  ||oo<  a0,  (2.2) 

where  is  a  constant  depending  on  the  flux. 

Recall  that  the  numerical  solutions  of  Godunov  scheme  are  the  cell  averages  of  the  exact 
solution  if  A  ||  (|w|  +  c)  ||oo<  1.  Thus  Lemma  2.2  implies  a0  =  1  for  the  Godunov  flux. 
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Following  the  same  proof  as  that  in  the  Appendix  of  [7],  it  is  straightforward  to  check  that 
tto  =  |  for  the  Lax-Friedrichs  flux. 


2.2  High  order  schemes 


We  now  consider  a  general  high  order  finite  volume  scheme,  or  the  scheme  satisfied  by  the 
cell  averages  of  a  DG  method  solving  (1.1),  with  Euler  forward  time  discretization,  which 
has  the  following  form 


wn+1  —  w”  \ 

wj  -  w)  “  A 


W 


■  1  ,  W  i 

3+2  J+2 


W  1  ,  W  .  1 

3~2 


(2.3) 


where  f  is  Godunov,  Lax-Friedrichs  or  kinetic  type  flux,  w "  is  the  approximation  to  the  cell 
average  of  the  exact  solution  v(x,t)  in  the  cell  Ij  =  [Xj_i,Xj+i]  at  time  level  n,  and  w7+1, 
are  the  high  order  approximations  of  the  point  values  v(Xj+i,tn)  within  the  cells  Ij  and 
Ij. |_i  respectively.  These  values  are  either  reconstructed  from  the  cell  averages  w"  in  a  hnite 
volume  method  or  read  directly  from  the  evolved  polynomials  in  a  DG  method.  We  assume 
that  there  is  a  polynomial  vector  qj(x)  =  (pj(x),rrij(x),  Ej(x))T  (either  reconstructed  in  a 
hnite  volume  method  or  evolved  in  a  DG  method)  with  degree  k,  where  k  >  1,  defined  on 
Ij  such  that  w”  is  the  cell  average  of  qj(x)  on  Ij,  =  qj(^-_i)  and  w7+i  =  clj(xj+ i)- 

We  need  the  Appoint  Legendre  Gauss-Lobatto  quadrature  rule  on  the  interval  Ij  = 
[xj_i,Xj+i\,  which  is  exact  for  the  integral  of  polynomials  of  degree  up  to  2N  —  3.  We  would 
need  to  choose  N  such  that  2N  —  3  >  k.  Denote  these  quadrature  points  on  Ij  as 


{x  i  =x1i,x2i,--  ,x 


:N—1 


3  ’  '  ~3  ’  X. 


■N 


=  X3+0- 


N 

Let  wa  be  the  quadrature  weights  for  the  interval  [— |,  |]  such  that  )T)  wa  —  1. 

a=  1 


(2.4) 


Theorem  2.3.  The  high  order  scheme  (2.3)  satisfies  a  minimum  entropy  principle,  i.e., 
assuming  the  numerical  solution  at  time  level  n  has  positive  density  and  positive  pressure, 
then  w ”+1  has  positive  density  and  positive  pressure,  and 


S( w"+1)  >  min  |min5'(qJ  ( 


<5? 


)).S(w++., 


J  2 
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under  the  CFL  condition 


A  II  (M  +  c)  ||oo<  Wia0. 

In  particular,  if  qj(xf)  G  G  for  all  j  and  a,  then  w ”+1  e  G. 


(2.5) 


Proof.  The  positivity  of  density  and  pressure  of  w"+  was  proved  in  Theorem  2.1  of  [13]. 


Thus  5'(w"+1)  is  well-defined.  The  exactness  of  the  quadrature  rule  for  polynomials  of  degree 


k  implies 


w  = 

3 


l  r  N 

—  /  qj(x)dx  =  '±2wa q 

■'  T:i  a= 1 


By  adding  and  subtracting  f  I  w+_i,w.  ,  ),  the  scheme  (2.3)  becomes 


N 


w?+1 


where 


«V>q.;  (•?,')  -  A 


a=l 


f(w7+i.w-+iJ-f(w-_i,w7+i 


f  W+  !,W.  ,  —  f  W.  1,W+  ! 

'  j-r  j+5/  V  j-r 


N—l 


J2™aq3' 


a=2 


A 

?j)  +  WN  1  W  1  -  — 
J  '  3  +  2  Wn 


3+2  3+1 


f  w  !,w;  !  -  f  wT,w  ! 


,  A 
+w i  (  w;  !  -  — 

3  2  Wi 

N-l 


f(w+  W.  ,  )  -  f  (W  1?W+  j 
'  3~r  3+^J  v  3~r  3-\ 


Y.  Wa<ij(x])  +  WNUN  +  WiHi, 


a=2 


3-r  j+i 


A 


Hi  =  w+  j  — 

3  2  Wi  L 


f(w+i,w.  f(w.  1.w+1 

1  3-h'  3+V  \  3~\'  3-\ 


Htv  =  w 


_A_ 

3+h  wn 


(2.6) 

(2.7) 


Notice  that  (2.6)  and  (2.7)  are  both  of  the  type  (2.1),  and  w\  =  wn,  therefore  Hi  and 
Htv  satisfy  the  minimum  entropy  principle  under  the  CFL  condition  (2.5).  Since  W“+1  is 
a  convex  combination  of  Hi,  Hjv  and  q by  Lemma  2.1,  we  get  the  minimum  entropy 


principle  for  w 


n+ 1 
3 


□ 


The  high  order  SSP  time  discretizations  will  keep  the  validity  of  Theorem  2.3  since  they 
are  convex  combinations  of  Euler  forward. 
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Theorem  2.3  implies  that,  to  have  the  minimum  principle  for  W"+1  G  G,  we  need  to 
enforce  q j(xf)  G  G.  The  positivity  of  density  and  pressure  of  q j(xf)  G  G  was  discussed  in 
[13].  Thus  here  we  only  show  how  to  enforce  the  entropy  part. 

At  time  level  n,  given  w"  G  G,  assume  q?(x“)  (a  —  1,  •  •  •  ,  N)  have  positive  density  and 
pressure.  Define  d G  =  {w  :  p,  p  >  0,  S  =  So  =  minx  5(w0(x))}  ,  and 


L(f)  =  (1  —  t)w'J  +  fqj(x),  0  <  t  <  1. 


(2.8) 


d G  is  a  surface  and  L(t)  is  the  line  segment  connecting  the  two  points  W"  and  q,,(a;),  where 
t  is  a  parameter.  If  S(qj(x))  <  So,  then  the  line  segment  L(f)  ( t  G  [0,1])  intersects  with 
the  surface  d G  at  one  and  only  one  point  since  G  is  a  convex  set.  If  S^q^a;))  <  So,  let  t(x) 
denote  the  parameter  in  (2.8)  corresponding  to  the  intersection  point;  otherwise  let  t(x )  =  1. 
In  practice,  we  can  find  t(x)  by  using  Newton  iteration  to  solve  S  (L(t(x)))  =  So-  Now  we 
define 


qj(z)  =  (q j(x)  -  w ”)  +  w",  0j  =  rnin  t{x).  (2.9) 

x£{xj,--- ,xf} 

The  limiter  (2.9)  should  be  used  for  each  stage  in  a  SSP  Runge-Kutta  method  or  each 
step  in  a  SSP  multi-step  method.  It  is  easy  to  check  that  the  cell  average  of  qj  (x)  over  is 


wj ' 


Lemma  2.4.  q,,(a;)  defined  in  (2.9)  satisfies  q jixf)  G  G  for  all  a. 


Proof.  First  notice  that  qj(a;)  is  a  convex  combination  of  qj(x)  and  w",  thus  q jixf)  still 
have  positive  density  and  pressure  since  pressure  is  a  concave  function.  We  only  need  to 
prove  S(qj(xf))  >  So  for  the  case  that  S(qj  (£“))  <  So¬ 
li  S(qLj(xf))  <  S0,  then  S  (L (t(xf)))  =  S0  and 


qj(^)  =  eJ  + 

e 


W; 


/(.?;)  (q'^)-  w^~w'-  V'  JTdf 


W; 


0i  -L«$?))+(l  6l 


w 


3  ’ 


t(xf)  V  t{xf 

So  qj  (xf)  is  a  convex  combination  of  L (t(xf))  and  wf,  thus  S(qj(xf))  >  So- 


□ 


We  refer  to  a  generic  smooth  solution  as  a  smooth  solution  v(x,fn )  G  G  satisfying 
min  S(v(x,  tn))  =  So  and  the  second  order  derivative  of  £(v(x,tn))  with  respect  to  x  does 
not  vanish  at  the  global  minimum.  For  such  generic  smooth  solutions,  the  limiter  (2.9) 
does  not  affect  the  high  order  accuracy  of  the  original  scheme.  Assume  q?(x)  is  a  (k  +  1)- 
th  order  accurate  approximation  qy  (x)  —  v(x,  tn)  =  0( Axk+1).  Without  loss  of  generality, 
assume  9j  =  t(xj)  for  some  (5.  Since  q L (t(xj))  and  w"  he  on  the  same  line,  we  have 


.  Thus, 


q j  (x)  -  qj  (x)  =  9j  (q,  (x)  -  w")  +  w?  -  q^  (x) 


=  (9j  -  1H<1; (.'•)- w") 
K'(S?)-L  (t(xf)) 


qj  (Xj )  -  w" 


W 


Define  d( z,  dG)  =  min  ||  z 

we9G 


w  || .  Since  So  is  the  minimum  of  S(v(x,tn)),  there 
is  at  least  one  x“  such  that  S0  —  S'(v(x“,  tn))  >  CAx2  where  C  is  a  nonzero  constant 


depending  on  the  derivatives  of  S(v(x,tn)).  This  implies  d(v(x",  tn),  dG)  >  0( Ax2),  thus 


N 


d(qj(x<j),dG)  >  0( Ax2),  w"  is  the  cell  average  of  q j(x)  implies  w”  =  waqj(x^).  So 


a=l 


d(Wj,  dG)  >  0( Ax2). 


On  the  other  hand,  Oj  =  t(Xj)  implies  S(qj(xj))  ^  G,  so  ||  q j(xj)  —  w"  ||>  d(wj,dG)  > 
0(Ax2).  The  overshoot  is  small  ||  q j(Xj)  —  L{t(xj))  ||=  0(Axk+1)  since  qj(x)  is  an  accurate 
approximation  to  v(x,  tn)  G  G. 

Finally,  notice  that  ||  qj(x)  —  w"  ||=  0( Ax),  we  get  that  ||  q j(x)  —  qj-(x)  ||=  0( Axk), 
Vx  G  Ij. 


We  remark  that  for  the  non-generic  situation  that  the  second  derivative  of  S'(v(x,  tn ))  with 
respect  to  x  does  vanish  at  the  global  minimum,  it  seems  difficult  to  design  a  conservative 
limiter  which  can  be  proved  not  to  destroy  accuracy.  On  the  other  hand,  the  fact  that  our 
limiter  is  easy  to  implement  also  for  multi- dimensional  cases  (see  next  section)  and  that  it 
maintains  high  order  accuracy  for  generic  smooth  solutions  makes  it  a  good  technique  to 
adopt  for  high  order  schemes. 
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3  The  two-dimensional  cases 


In  this  section  we  extend  our  result  to  finite  volume  or  DG  schemes  of  (k  +  l)-th  order 
accuracy  solving  two-dimensional  Euler  equations  with  initial  data  wq (x,y) 


+  f(w)x  +  g(w)y  =  0,  t  >  0,  (x,y)  E  R2, 


,  f(w)  = 


pu2  +  p 
puv 

(E  +  p)u 


g(w )  = 


pir  +  p 
(E  +  p)v 


where  m  =  pu,  n  =  pv,  E  =  ^ pu 2  +  \pv2  +  pe,p  =  (7  —  1  )pe,  and  ( u ,  v)  is  the  velocity.  The 

eigenvalues  of  the  Jacobian  f'(w)  are  u  —  c,u,u  and  u  +  c  and  the  eigenvalues  of  the  Jacobian 

g'(w)  are  v  —  c,  v,  v  and  v  +  c.  The  specific  entropy  S  =  In  is  quasi-concave  with  respect  to 

w  if  p  >  0  and  the  set  of  admissible  states  G  =  <  w|  p  >  0,p  >  0,  S  >  So  =  min  i5'(w0(x,  y ))  > 

[  x,y  J 

is  still  convex. 


3.1  Rectangular  meshes 

For  simplicity  we  assume  we  have  a  uniform  rectangular  mesh.  At  time  level  n,  we  have  an 
approximation  polynomial  q ij(x,y)  of  degree  k  with  the  cell  average  w”.  on  the  ( i,j )  cell 
[Xi_  i,xi+i]  x  [y^yj+i].  Let  w+  ^(y),  w~+  ^.(y),  w+_i(x),  w“+i(a:)  denote  the  traces  of 
q ij(x,y)  on  the  four  edges  respectively.  A  finite  volume  scheme  or  the  scheme  satisfied  by 
the  cell  averages  of  a  DG  method  on  a  rectangular  mesh  can  be  written  as 


where  f(-,  •),  g(-,  •)  are  one  dimensional  fluxes.  The  integrals  can  be  approximated  by  quadra¬ 
tures  with  sufficient  accuracy.  Let  us  assume  that  we  use  a  Gauss  quadrature  with  L  points, 
which  is  exact  for  single  variable  polynomials  of  degree  k.  We  assume  Sf  =  {xf  :  f3  = 
1,  •  •  •  ,  L)  denote  the  Gauss  quadrature  points  on  [xj_i ,  xi+i],  and  AJ  =  {y^  :  (3  —  1,  •  •  •  ,  L} 
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denote  the  Gauss  quadrature  points  on  [yj_i,yj+i\.  For  instance,  (aq_i  ,y?)  {(3  —  1,  •  •  •  ,  L) 
are  the  Gauss  quadrature  points  on  the  left  edge  of  the  (i,j)  cell.  The  subscript  (3  will 
denote  the  values  at  the  Gauss  quadrature  points,  for  instance,  ^  =  ~wt_i  -{Vj)-  Also, 
W/3  denotes  the  corresponding  quadrature  weight  on  interval  [— |,  |],  so  that  Y^3=iwP  =  1- 
We  will  still  need  to  use  the  iV-point  Gauss-Lobatto  quadrature  rule  where  N  is  the  smallest 
integer  satisfying  2N  —  3  >  k,  and  we  distinguish  the  two  quadrature  rules  by  adding  hats 
to  the  Gauss-Lobatto  points,  i.e.,  Sf  =  {xf  :  a  —  1,  •  •  •  ,  N}  will  denote  the  Gauss-Lobatto 
quadrature  points  on  [,x'8_i ,  xi+3\,  and  Sy  =  { y°f  :  a  —  1,  •  •  •  ,N}  will  denote  the  Gauss- 
Lobatto  quadrature  points  on  [yj_i,Uj+i].  Subscripts  or  superscripts  f3  will  be  used  only  for 
Gauss  quadrature  points  and  a  only  for  Gauss-Lobatto  points. 

Let  A[  =  Ai  and  \2  =  AL_  then  the  scheme  becomes 

L 


w 


■71+1 

ij 


=  W 


0=1 


fK+i 


^wi+b0> 


f(w  ■  i 


,0^U,0> 


Y 

0=1 


g(W 


00+V 


w 


00+h' 


g(w 


00-V 


w 


0j-\‘ 


(3.2) 


We  use  <g)  to  denote  the  tensor  product,  for  instance,  Sf  ®  Sy  =  {(x,  y)  :  x  G  Sf,  y  G  Sy}. 
Define  the  set  Sij  as 

Sii  =  (Sf®S>)  U  (3.3) 


For  simplicity,  let  /i  =  arid  /w  =  where  a,  =||  (|u|  +  c)  a2  =| 

(M  +  c)  || oo-  Notice  that  Wi  =  wn,  we  have 


_ n  Ml 

w  •  =  - 

13  AxAy 


rxi+i  fyi+l 


q  ij(x,y)dydx  + 


M2 


ryi+h  fx*+h 


’  x.  i  J y .  i 


y .  i  ox.  i 


v?) 


AxAy 

L  N  L  N 

=  l-l  \  EE  wpwaq  ij(xf,yj)  +  M2  EE  wpwaqij(x f, 

/3=1  o=l  /3=1  o=l 

L  IV— 1 

=  EE  |//iqy(x?,j/f)  +  /biqp(a+ ,?/ 

0=1  o=2 
L 

+ Y 

0=1 


q  ij(x,y)dxdy 


3°) 

j  > 


+  fiw+i  „  +  .+,  +  p->wy_. 


(3.4) 
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Theorem  3.1.  Consider  a  two-dimensional  finite  volume  scheme  or  the  scheme  satisfied 
by  the  cell  averages  of  the  DG  method  on  rectangular  meshes  (3.2),  associated  with  the 
approximation  polynomials  q ij(x,y)  of  degree  k  (either  reconstruction  or  DG  polynomials) . 
If  wtj±1’  wi±i  g  e  ^  and  £  G  (for  any  (x,y)  G  Sij),  then  w"+1  G  G  under  the 

CFL  condition 

Aiai  T  A2O2  —  ^ '  1  o 0 ■  (3.5) 


Proof.  Plugging  (3.4)  in,  (3.2)  can  be  written  as 

L  N-l 

fi  _  ' 

’  o 


W”+1  =  22^2  W^° 


0=1  a= 2 
L 

+1*1  22  UJfjVJi 
0=1 


/pq Vj)  +  y) 

Ai 


w 


*+P/3 


PlWi 


f(wi+i  p,  W+  1  p  -  f(w+  1  „,  w-  i 


i-h0’  i+G0' 


w 


Ai 


Tw^'wAd~f(wA?'w.-i.. 


T  1*2  22  W0™N 
0=1 


A2 


W 


Ai+I 


/*2'UJ\ 


^w0,Mfiw0oH>-^w 


0d-Gw0,j+G 


+  A2 

1-  — ieiw 


/D-P  W/hi+D 


"2  Ai2Wl 

Following  the  same  arguments  as  in  Theorem  2.3,  we  conclude  w”+1  G  G. 


□ 


The  limiter  in  the  previous  section  can  be  extended  easily  to  two-dimensional  cases.  At 
time  level  n,  given  w”  G  G,  do  the  following  modification 


q ij(x,  y)  =  Oij  (q ij(x,  y)  -  W T)  +  w"-,  %  =  min  t(x,  y), 

( x,y)eSij 


(3.6) 


where  t(x,  y)  is  the  parameter  corresponding  to  the  intersection  point  of  the  surface  d G  and 
the  line  segment  L(f)  =  (1  —  f)w™-  +  tqtfix,  y)  if  q?J(x,  y)  G ;  t(x,y )  =  1  otherwise. 


3.2  Triangular  meshes 


For  simplicity,  we  only  discuss  DG  schemes  in  this  subsection.  All  the  conclusions  will  also 
hold  for  a  high  order  finite  volume  scheme. 
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For  each  triangle  K  we  denote  by  VK  (i  =  1,  2,  3)  the  length  of  its  three  edges  elK  (?'  = 
1,2,3),  with  outward  unit  normal  vector  vl  ( i  =  1,2,3).  K{i)  denotes  the  neighboring 
triangle  along  4  and  \K\  is  the  area  of  the  triangle  K .  Let  F(w,  v,  v)  be  a  one  dimensional 
monotone  flux  in  the  v  direction  satisfying  F(w,v,  v)  =  —  F(v,  w,  —  u)  (conservation),  and 
F(w,  w,  v)  =  F(w)  •  v  (consistency),  with  F(w)  =  (f(w),g(w)).  For  example,  the  Lax- 
Friedrichs  flux  is  dehned  as 


F(w,  v,  v)  =  -(F(w)  ■  i/  +  F(v)  •  v-  a(v  -  w)),  a  =||  \{u,v)\  +  c 


A  high  order  finite  volume  scheme  or  a  scheme  satisfied  by  the  cell  averages  of  a  DG 
method,  with  first  order  Euler  forward  time  discretization,  can  be  written  as 


TTT«  +  1  _ 

TK 


At 


=W 


F(w 


int.(K)  ext(K) 

.  W, 


,  z/*)ds, 


i=  1  ^ K 


where  is  the  cell  average  over  K  of  the  numerical  solution,  and  'w>)nt<K\  w('  xt('K'>  are  the 
approximations  to  the  values  on  the  edge  eK  obtained  from  the  interior  and  the  exterior  of 
K.  Assume  the  DG  polynomial  on  the  triangle  K  is  q/v(rc,  y)  of  degree  k ,  then  in  the  DG 
method,  the  edge  integral  should  be  approximated  by  the  ( k  +  l)-point  Gauss  quadrature. 
The  scheme  becomes 


w 


n+ 1 
K 


=  W 


K 


At  V- 

1  1  i= 1  /3=  1 


w 


int.(K)  ext(K) 


i,/3 


W 


i,f3 


\ul)w0l 


i 

K  i 


(3.7) 


where  Wp  denote  the  ( k  +  l)-point  Gauss  quadrature  weights  on  the  interval  [—A,  A],  so 

k-\-l 

that  Y2  w/3  —  1,  and  and  wf-xl,(K}  denote  the  values  of  w  evaluated  at  the  /1-th 

0=  i 


Gauss  quadrature  point  on  the  i-th  edge  from  the  interior  and  exterior  of  the  element  K 
respectively. 

We  need  the  quadrature  rule  introduced  in  [14]  for  q k(x,i/)  on  K .  In  the  barycentric 
coordinates,  the  set  of  quadrature  points  for  polynomials  of  degree  k  on  a  triangle  K  can 
be  written  as 


SkK 


+  £°)( 


1 

2 


5 
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(3.8) 


((j -«P)(5-A5  +  «'.(5  +  8*)(5 -«')). 

((i  +  Sf)(5-A  + 

where  ua  (a  =  1,  •  •  •  ,N)  and  A9  (/3  =  l,---  ,k  +  1)  are  the  Gauss-Lobatto  and  Gauss 
quadrature  points  on  the  interval  [— |,  |]  respectively.  See  Figure  3.1  for  an  illustration  of 

si 


Figure  3.1:  Points  in  (3.8)  for  k  =  2. 

Following  Theorem  5.1  in  [14]  and  the  same  arguments  as  in  Theorem  2.3,  we  have 

Theorem  3.2.  For  the  scheme  (3.7)  with  the  polynomial  qxix,  y )  of  degree  k,  if  G  G 

and  qK(x,y)  G  G,  \/(x,y)  G  S ^  where  S'f  is  defined  in  (3.8),  then  w^G1  g  G  under  the  CFL 

3 

condition  aM |  Ir  —  fwicco- 

i= 1 

The  the  same  limiter  can  be  used  to  enforce  the  sufficient  condition.  At  time  level  n, 
given  w \  G  G,  do  the  following  modification 

q k(x,  y)  =  0K  (qK(x,  y)  -  w£)  +  w£,  0K  =  min  t(x,  y),  (3.9) 

( x,y)eSK 

where  t(x,  y)  is  the  parameter  corresponding  to  the  intersection  point  of  the  surface  d G  and 
the  line  segment  L(f)  =  (1  —  t)WlK  +  tqK(x,y)  if  q K(x,y)  ^  G ;  t(x,y )  =  1  otherwise. 

4  Numerical  Tests 

Example  4.1.  Accuracy  tests. 
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We  first  test  the  accuracy  of  the  limiter  (2.9).  The  initial  condition  is  p0(x,y)  =  1  + 
|sin(27nr),  Mo(^)  =  1,  Po(x)  =  1.  The  domain  is  [0,1]  and  the  boundary  condition  is  peri¬ 
odic.  The  exact  solution  is  p(x,y,t )  =  1  +  |sin(27r(a;  —  t )),  u(x,t )  =  1,  p(x,t)  =  1.  The 
numerical  schemes  are  the  third  order  and  fourth  order  DG  schemes  with  Lax- Friedrichs  flux 
[2]  and  the  third  order  SSP  time  discretizations  with  the  limiter  (2.9)  used  at  each  time  stage 
or  each  time  step. 

The  third  order  SSP  Runge-Kutta  method  in  [9]  (with  the  CFL  coefficient  c  =  1)  is 

«(1)  =  un+AtF(un) 

w(2)  =  -«"  +  -(««  +  AtF(uW) 

4  4 

un+l  =  i„"  A/F((/(2))) 

o  o 

where  F(u)  is  the  spatial  operator,  and  the  third  order  SSP  multi-step  method  in  [8]  (with 
the  CFL  coefficient  c  =  |)  is 

f”+1  =  +  3A«F(u”))  +  K"3)). 

Here,  the  CFL  coefficient  c  for  a  SSP  time  discretization  refers  to  the  fact  that,  if  we  assume 

the  Euler  forward  time  discretization  for  solving  the  equation  Ut  =  F(u )  is  stable  in  a  norm 

or  a  semi-norm  under  a  time  step  restriction  At  <  A t0l  then  the  high  order  SSP  time 
discretization  is  also  stable  in  the  same  norm  or  semi-norm  under  the  time  step  restriction 
At  <  cAt0. 

For  k  —  2  and  k  —  3,  wi  —  |  and  a0  =  thus  the  time  step  (2.5)  for  P 2  DG  with 
Runge-Kutta  is  taken  as  At  =  ,m  .^L..  .  Since  the  CFL  coefficient  c  —  \  for  the  third 
order  SSP  multi-step  method,  the  time  step  is  taken  as  At  =  .  For  the  fourth 

order  scheme,  in  order  to  make  the  error  from  spatial  discretizations  dominant,  we  replace 
Ax  by  Aa;3. 

The  accuracy  result  is  listed  in  Table  4.1.  We  observe  that,  for  Runge-Kutta,  the  accu¬ 
racy  degenerates  when  the  mesh  is  fine  enough.  This  is  due  to  the  lower  order  accuracy  in 
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the  intermediate  stages  of  the  Runge-Kutta  method.  In  particular,  recall  that  the  limiter 
(2.9)  does  not  destroy  accuracy  for  generic  smooth  solutions  only  if  the  polynomial  qj(x) 
is  a  (k  +  l)-th  accurate  approximation  to  the  exact  solution.  The  DG  polynomials  qj(x) 
in  the  intermediate  stages  of  a  Runge-Kutta  methods  are  in  general  not  ( k  +  l)-th  order 
accurate,  therefore  the  limiter  (2.9)  may  kill  the  accuracy  when  it  is  imposed  in  the  interme¬ 
diate  stages.  The  same  phenomenon  exists  for  the  high  order  maximum-principle-satisfying 
schemes,  see  [12],  A  similar  phenomenon  of  the  Runge-Kutta  method  in  the  context  of 
boundary  conditions  was  pointed  out  in  [1], 

The  full  accuracy  order  is  observed  for  the  multi-step  time  discretization,  which  justifies 
that  the  limiter  itself  does  not  kill  accuracy.  Since  accuracy  degeneracy  is  usually  only 
observed  on  very  fine  meshes  for  Runge-Kutta  methods,  in  applications  it  is  often  acceptable 
to  use  the  Runge-Kutta  methods,  similar  to  the  conclusions  in  [1,  12], 


Table  4.1:  Third  order  SSP  time  discretizations  and  high  order  DG  spatial  discretizations 
with  the  limiter  (2.9),  Ax  =  A,  t=0.1. 


N 

SSP  Runge-Kutta 

SSP  multi-step 

L1  error 

order 

L°o  error 

order 

L1  error 

order 

L°°  error 

order 

P2DG 

At  = 

1.49E-3 

l 

Ax 

At  = 

1.58E-3 

i 

Ax 

8 

12  ||(|u 

l+c)||oo 

5.23E-3 

_ 

36  ||(|ii 

I+C)||oo 

5.79E-3 

_ 

16 

1.64E-4 

3.17 

8.46E-4 

2.62 

1.65E-4 

3.25 

8.77E-4 

2.72 

32 

2.07E-5 

2.98 

9.07E-5 

3.22 

2.06E-5 

3.00 

9.07E-5 

3.27 

64 

2.62E-6 

2.97 

1.17E-5 

2.95 

2.60E-6 

2.98 

1.17E-5 

2.95 

128 

3.30E-7 

2.98 

1.47E-6 

2.98 

3.25E-7 

2.99 

1.47E-6 

2.98 

256 

4.16E-8 

2.99 

1.84E-7 

2.99 

4.07E-8 

2.99 

1.84E-7 

2.99 

512 

5.23E-9 

2.99 

3.51E-8 

2.39 

5.09E-9 

3.00 

2.31E-8 

3.00 

1024 

6.60E-10 

2.98 

1.02E-8 

1.78 

6.36E-10 

3.00 

2.88E-9 

3.00 

P3DG 

At  = 

5.11E-4 

3 

1  Ax? 

At  = 

5.19E-4 

3 

1  Ax? 

4 

12  ||(|u 

l+«0||oo 

4.52E-3 

_ 

36  ||(|u 

l+c)||oo 

2.09E-3 

_ 

8 

2.45E-5 

4.38 

1.05E-4 

4.31 

2.46E-5 

4.39 

1.05E-4 

4.31 

16 

1.40E-6 

4.12 

5.82E-6 

4.18 

1.40E-6 

4.13 

5.34E-6 

4.30 

32 

9.02E-8 

3.96 

5.61E-7 

3.38 

8.41E-8 

4.06 

3.74E-7 

3.83 

64 

6.66E-9 

3.75 

1.24E-7 

2.16 

5.21E-9 

4.01 

2.23E-8 

4.06 

128 

5.29E-10 

3.65 

1.83E-8 

2.77 

3.24E-10 

4.00 

1.42E-9 

3.97 

256 

4.37E-11 

3.59 

3.39E-9 

2.43 

2.02E-11 

4.00 

9.00E-11 

3.98 

512 

4.14E-12 

3.39 

5.61E-10 

2.59 

1.26E-12 

3.99 

5.55E-12 

4.01 

16 
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(d)  With  (2.9) 


X 
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X 

(f)  With  (2.9) 


Figure  4.1:  Lax  shock  tube  problem.  P3  element  DG  with  the  positivity-preserving  limiter 
without  any  TVD  or  TVB  limiter.  The  solid  lines  are  the  exact  solutions.  The  symbols  are 
the  numerical  solutions.  ^ 


Example  4.2.  Lax  shock  tube  problem. 


For  high  order  DG  schemes  solving  the  compressible  Euler  equations,  even  though  the 
characteristicwise  TVB  limiter  in  [2]  can  kill  oscillations,  it  is  not  sufficient  to  stabilize 
the  scheme  for  problems  with  low  densities  or  low  pressures.  In  [13],  it  was  reported  that 
high  order  RKDG  scheme  with  both  the  positivity-preserving  limiter  and  the  TVB  limiter 
worked  fine  for  a  lot  of  demanding  problems.  Recent  study  reveals  that  the  third  order 
RKDG  scheme  with  only  the  positivity-preserving  limiter  is  stable  even  for  strong  shocks, 
see  [11],  which  is  not  surprising  since  a  conservative  positivity-preserving  scheme  is  L 1  stable 
[15].  However,  the  positivity-preserving  limiter  alone  can  not  kill  oscillations  for  high  order 
DG  schemes,  and  the  oscillations  are  much  more  prominent  in  the  fourth  order  DG  scheme 
than  in  the  third  order  scheme.  See  Figure  4.1  for  the  result  of  the  fourth  order  DG  scheme 
with  P3  element  for  the  Lax  shock  tube  problem.  We  can  see  that  the  result  with  only 
the  positivity-preserving  limiter  is  oscillatory  and  the  result  with  the  positivity-preserving 
limiter  and  (2.9)  has  much  less  oscillations.  In  other  words,  enforcing  the  minimum  entropy 
principle  will  damp  the  oscillations  in  high  order  schemes,  which  was  pointed  out  in  [6]. 

5  Concluding  Remarks 

In  this  paper,  we  have  discussed  the  minimum  entropy  principle  for  high  order  schemes  solv¬ 
ing  the  compressible  Euler  equations  in  gas  dynamics.  An  extension  of  positivity-preserving 
limiter  in  [13]  can  be  used  to  enforce  the  minimum  entropy  principle.  The  generalizations 
to  higher  dimension  are  straightforward.  Numerical  tests  imply  that  enforcing  the  minimum 
entropy  principle  may  damp  the  oscillations  in  high  order  schemes. 

Acknowledgments:  The  authors  would  like  to  thank  Eitan  Tadmor  for  helpful  discussions 
about  the  minimum  entropy  principle  for  gas  dynamics  equations. 
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